1 데이터와 예측모형12

앞서 작업한 데이터와 선형회귀모형과 ARIMA 모형을 저장한 data/corona_fcst_list.rds 리스트 객체를 단변량 시계열 분석 작업을 위해서 풀어둔다.

library(tidyverse)
library(tidymodels)

library(timetk)
library(modeltime)

# 데이터 + 모형 ----
corona_fcst_list <- read_rds("data/corona_fcst_list.rds")

# 데이터 ----
full_tbl <- corona_fcst_list$data

# 모형 ----
wkfl_fit_lm <- corona_fcst_list$model$wkfl_fit_lm
wkfl_fit_arima <- corona_fcst_list$model$wkfl_fit_arima

1.1 훈련/시험 데이터 분할

예측할 데이터와 모형개발에 활용할 데이터로 나눈 후에 모형개발에 활용할 데이터를 훈련/시험 데이터로 나눈다.

## 예측 데이터와 모형개발 데이터로 분리
forecast_tbl <- full_tbl %>% 
  filter(is.na(확진자))

history_tbl <- full_tbl %>% 
  filter(!is.na(확진자))

## 훈련/시험 데이터 분할

splits <- history_tbl %>% 
  time_series_split(date_var    = 날짜,  
                    assess      = 30,
                    cumulative  = TRUE)
splits
<Analysis/Assess/Total>
<309/30/339>

훈련/시험 데이터로 잘 나눠졌는지 시각적으로 확인한다.

splits %>% 
  tk_time_series_cv_plan() %>% 
  plot_time_series_cv_plan(.date_var = 날짜,
                            .value   = 확진자)

2 피처 공학

tidymodels 생태계의 recipe 팩키지 recipe() 함수를 사용해서 피쳐 공학(feature engineering) 작업을 수행한다. 추후, 다양한 피처를 작업할 예정이라 대략적인 틀만 잡아둔다.

recipe_spec <- recipes::recipe(확진자 ~ 날짜, data = training(splits))

recipe_spec %>% prep() %>%  juice()
# A tibble: 309 x 2
   날짜       확진자
   <date>      <int>
 1 2020-01-23      0
 2 2020-01-24      1
 3 2020-01-25      0
 4 2020-01-26      1
 5 2020-01-27      1
 6 2020-01-28      0
 7 2020-01-29      0
 8 2020-01-30      0
 9 2020-01-31      7
10 2020-02-01      1
# ... with 299 more rows

3 단변량 예측 알고리즘

단변량 예측 알고리즘으로 다음 세가지 대표적인 알고리즘을 적용시킨다.

  • ARIMA
  • ETS
  • TBATS
  • Prophet

3.1 ARIMA

대표적인 시계열 데이터 예측모형 ARIMA로 적합시킨다.

mdl_arima_spec <- arima_reg() %>%
    set_engine("auto_arima")

wkfl_fit_arima <-  workflows::workflow() %>% 
  add_recipe(recipe_spec) %>% 
  add_model(mdl_arima_spec) %>% 
  fit(training(splits))

wkfl_fit_arima
== Workflow [trained] ==========================================================
Preprocessor: Recipe
Model: arima_reg()

-- Preprocessor ----------------------------------------------------------------
0 Recipe Steps

-- Model -----------------------------------------------------------------------
Series: outcome 
ARIMA(0,1,1) 

Coefficients:
          ma1
      -0.2198
s.e.   0.0519

sigma^2 estimated as 2729:  log likelihood=-1654.94
AIC=3313.87   AICc=3313.91   BIC=3321.33

3.2 ETS

지수평활(exponential smoothing) 중 ETS(Error-Trend-Season)를 예측모형으로 적합시킨다.

mdl_ets_spec <- exp_smoothing(
    error   = "auto",
    trend   = "auto",
    season  = "auto",
    damping = "auto"
) %>%
    set_engine("ets")

wkfl_fit_ets <-  workflows::workflow() %>% 
  add_recipe(recipe_spec) %>% 
  add_model(mdl_ets_spec) %>% 
  fit(training(splits))

wkfl_fit_ets
== Workflow [trained] ==========================================================
Preprocessor: Recipe
Model: exp_smoothing()

-- Preprocessor ----------------------------------------------------------------
0 Recipe Steps

-- Model -----------------------------------------------------------------------
ETS(A,N,A) 

Call:
 forecast::ets(y = outcome, model = model_ets, damped = damping_ets) 

  Smoothing parameters:
    alpha = 0.7826 
    gamma = 1e-04 

  Initial states:
    l = 2.4739 
    s = 1.6161 1.1827 -17.5283 -11.5425 14.7329 4.8602
           6.6788

  sigma:  51.3469

     AIC     AICc      BIC 
4216.526 4217.265 4253.860 

3.3 TBATS

TBATS를 예측모형으로 적합시킨다.

mdl_tbats_spec <- seasonal_reg(
  mode              = "regression",
  seasonal_period_1 = "auto"
) %>%
    set_engine("tbats")

wkfl_fit_tbats <-  workflows::workflow() %>% 
  add_recipe(recipe_spec) %>% 
  add_model(mdl_tbats_spec) %>% 
  fit(training(splits))

wkfl_fit_tbats
== Workflow [trained] ==========================================================
Preprocessor: Recipe
Model: seasonal_reg()

-- Preprocessor ----------------------------------------------------------------
0 Recipe Steps

-- Model -----------------------------------------------------------------------
TBATS(1, {0,0}, -, {<7,2>})

Call: forecast::tbats(y = outcome, use.parallel = use.parallel)

Parameters
  Alpha: 0.7812103
  Gamma-1 Values: -0.002287667
  Gamma-2 Values: 0.00209979

Seed States:
          [,1]
[1,] 21.234291
[2,] 11.179662
[3,] -6.490442
[4,]  2.899411
[5,] -3.565564

Sigma: 50.80503
AIC: 4215.104

3.4 Prophet

페이스북에서 개발한 Prophe 예측모형으로 적합시킨다.

mdl_prophet_spec <- prophet_reg(
) %>%
    set_engine("prophet")

wkfl_fit_prophet <-  workflows::workflow() %>% 
  add_recipe(recipe_spec) %>% 
  add_model(mdl_prophet_spec) %>% 
  fit(training(splits))

wkfl_fit_prophet
== Workflow [trained] ==========================================================
Preprocessor: Recipe
Model: prophet_reg()

-- Preprocessor ----------------------------------------------------------------
0 Recipe Steps

-- Model -----------------------------------------------------------------------
PROPHET Model
- growth: 'linear'
- n.changepoints: 25
- changepoint.range: 0.8
- yearly.seasonality: 'auto'
- weekly.seasonality: 'auto'
- daily.seasonality: 'auto'
- seasonality.mode: 'additive'
- changepoint.prior.scale: 0.05
- seasonality.prior.scale: 10
- holidays.prior.scale: 10
- logistic_cap: NULL
- logistic_floor: NULL
- extra_regressors: 0

4 모형 평가

workflow 객체를 modeltime 팩키지 modeltime_table() 함수에 넣어 객체로 만든 후에, 추가로 지수평활법 EST 예측모형을 추가한다. modeltime_accuracy() 함수를 사용해서 모형 성능을 파악한다.

model_tbl <- modeltime_table(
    wkfl_fit_arima,
    wkfl_fit_ets,
    wkfl_fit_tbats,
    wkfl_fit_prophet
    ) %>% 
  update_model_description(.model_id = 1, "ARIMA") %>% 
  update_model_description(.model_id = 2, "ETS") %>% 
  update_model_description(.model_id = 3, "TBATS") %>% 
  update_model_description(.model_id = 4, "Prophet")

model_tbl %>% 
  modeltime::modeltime_accuracy(testing(splits))
# A tibble: 4 x 9
  .model_id .model_desc .type   mae  mape  mase smape  rmse     rsq
      <int> <chr>       <chr> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>
1         1 ARIMA       Test   275.  30.3  2.78  36.7  341. NA     
2         2 ETS         Test   277.  30.5  2.80  37.0  343.  0.0720
3         3 TBATS       Test   277.  30.5  2.81  37.1  343.  0.0650
4         4 Prophet     Test   603.  73.0  6.10 116.   647.  0.461 

4.1 시각화

시험데이터를 통해 모형으로 예측한 값을 시각화한다. ETS MAE 값이 ARIMA 예측모형과 별다른 차이가 나고 있지 않지만, 신뢰구간이 생긴 것은 그나마 다행이다.

calibration_tbl <- model_tbl %>% 
  modeltime_calibrate(
    new_data = testing(splits)
  )

calibration_tbl %>% 
    modeltime_forecast(
        new_data      = testing(splits),
        actual_data   = history_tbl,
        keep_data     = TRUE,
        conf_interval = 0.10
    ) %>%
    plot_modeltime_forecast(
        .legend_max_width    = 60,
        .legend_show         = TRUE,
        .conf_interval_show  = TRUE,
        .conf_interval_alpha = 0.20,
        .conf_interval_fill  = "lightblue",
        .title = "코로나19 확진자 1개월 예측"
    )

5 확진자 예측

앞선 모형은 history_tbl 을 훈련/시험 데이터에 대해 적합을 시킨 것이라 … 이제 시간을 확대하여 modeltime_refit() 함수를 사용해서 모형을 전체 데이터에 대해 다시 적합시킨다.

refit_tbl <- calibration_tbl %>% 
  modeltime_refit(data = history_tbl)

마지막으로 앞서 구축된 모형을 바탕으로 현재까지 입수된 데이터를 바탕으로 한달 후를 예측한다.

refit_tbl %>% 
  modeltime_forecast(
    new_data    = forecast_tbl,
    actual_data = history_tbl,
    conf_interval = 0.3
  )  %>%
    plot_modeltime_forecast(
        .legend_max_width = 25,
        .conf_interval_fill = "lightblue",
        .conf_interval_alpha = 0.7,
        .interactive = TRUE
    )

6 모형과 데이터 저장

마지막 단계로 데이터와 모형을 저장하여 다음 단계로 나가기 위한 작업을 수행한다.

corona_fcst_list <- list(
  
  # 데이터 ---
  data = full_tbl, 
  
  # 예측모형 ----
  model = list(
    wkfl_fit_lm      = wkfl_fit_lm,
    wkfl_fit_arima   = wkfl_fit_arima,
    wkfl_fit_ets     = wkfl_fit_ets,
    wkfl_fit_tbats   = wkfl_fit_tbats,
    wkfl_fit_prophet = wkfl_fit_prophet
  )
)

# 모형 + 데이터 저장 ----
corona_fcst_list %>% 
  write_rds("data/corona_fcst_list.rds")
 

데이터 과학자 이광춘 저작

kwangchun.lee.7@gmail.com